Exact diagonalization results for an anharmonically trapped Bose-Einstein condensate 
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We consider bosonic atoms that rotate in an anharmonic trapping potential. Using numerical 
diagonalization of the Hamiltonian, we identify the various phases of the gas as the rotational 
frequency of the trap and the coupling between the atoms are varied. 
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I. INTRODUCTION 



The behavior of confined, Bose-Einstein condensed 
atoms under rotation has been studied extensively in re- 
cent years, both experimentally 0, 0, H El HI j as weu as 
theoretically. One factor which is very crucial in these 
studies is the form of the confinement. In most of the ex- 
periments performed on cold gases of atoms, the trapping 
potential is harmonic. It is interesting, however, that 
the harmonic trapping is special in many respects, as we 
show below. For this reason recent theoretical studies 

i0Biiiaiiiiiiiiiiiiiii|iEiiiiaiiai2ii, as 

well as the experiment of Ref. [231 > have considered traps 
with other functional forms, in which the trapping poten- 
tial grows more rapidly than quadratically at distances 
far away from the center of the cloud. Such trapping 
potentials introduce many novel phases. 

Motivated by this observation, we examine the low- 
est state of a Bose-Einstein condensate in a rotating an- 
harmonic trap. Up to now, all theoretical studies that 
have been performed on this problem were based upon 
the mean-field approximation. However, especially in the 
case of effective attractive interactions between the atoms 
[23L |24| , the use of the mean-field approximation is ques- 
tionable, as we explain in more detail below. In this 
study, we present "exact" results - for small numbers of 
atoms - from numerical diagonalization of the Hamilto- 
nian, which goes beyond the mean-field approximation. 

Experimentally it is possible to tune both the fre- 
quency of rotation of the trap, as well as the strength 
of the interatomic interaction. Motivated by these facts, 
we consider here a fixed trapping potential and examine 
the various phases and the corresponding phase diagram 
as function of the rotational frequency and of the cou- 
pling. 

In what follows we present our model in Sec. II. In 
Sec. Ill we discuss the general structure of the phase dia- 
gram. In Sec. IV we describe the details of the numerical 
diagonalization. In Sees. V and VI we present and ana- 
lyze our results for effective repulsive and attractive in- 
teratomic interactions, respectively. Finally, in Sec. VII 
we summarize our main results. 



II. MODEL 

The Hamiltonian that we consider consists of the 
single-particle part /i(rj), and of the interaction V( r i — 
Tj), which is taken to be of the form of a contact poten- 
tial, V(y l - Tj) = U 6(ri - r,-), 
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Here N is the number of atoms, and Uq = ^itTi 2 a/M, 
where a is the scattering length for elastic atom-atom 
collisions, and M is the atom mass. The single-particle 
part of the Hamiltonian has the usual form 



Mr) = -^rjV 2 + V tI . ap (r), 



(2) 



where the trap is assumed to be symmetric around the 
z axis. This is also taken to be the axis of rotation. In 
cylindrical coordinates p, <f>, and z, 
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Here uj is the trap frequency of the harmonic part of 
Vtrap, A is a dimensionless constant which is the coeffi- 
cient of the anharmonic term in Vt rap (taken to be much 
smaller than unity in this study, as is also the case in the 
experiment of Ref. 22] ) , and ao = y/h/Mco is the oscilla- 
tor length. Along the z axis we assume that the trapping 
potential V(z) is sufficiently strong that the typical sepa- 
ration between single-particle energy levels is much larger 
than the typical interaction energy. This assumption im- 
plies that the motion along the z axis is frozen out and 
our problem is essentially two-dimensional. 



III. GENERAL STRUCTURE OF THE PHASE 
DIAGRAM 

As mentioned in the introduction, motivated by the 
experimental situation, we here examine the phase di- 
agram as function of the rotational frequency f2 of the 
trap and of the strength of the interaction. It is natural 
to measure f2 in units of u and the strength of the inter- 
action ~ nUo, where n is the atomic density, in units of 
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the oscillator energy hoj. The corresponding dimension- 
less ratio nUo/fiLu is ~ era, where a = N/Z is the density 
per unit length. Here Z is the width of the atomic cloud 
along the z axis. 

It is instructive to first consider the general structure 
of the phase diagram 0, 0, ^| . This consists of several 
distinct phases. Let us start with repulsive interactions, 
a > 0. For a fixed il, as era increases, the total energy 
eventually becomes dominated by the interaction energy 
which is minimized via the formation of singly quantized 
vortex states. This is similar to vortex formation in liquid 
Helium, or in harmonically trapped gases. 

On the other hand, for a constant interaction strength, 
as the rotational frequency increases - and the single- 
particle part is the dominant part of the Hamiltonian - 
the total energy of the gas is minimized via formation 
of multiply quantized vortices. This is most easily seen 
if one thinks of the effective potential felt by the atoms, 
which consists of Vtrap minus the centrifugal potential, 



Kff(r) 



^trap(r) - -MQ 2 p 2 = 

= ^ 2 -flV + ^^ + y(4 (4) 

For the form of V tTap we have considered, V e g takes the 
form of a mexican-hat shape for Q > w, and thus the 
atoms prefer to reside along the bottom of this poten- 
tial. This is precisely the density distribution that corre- 
sponds to multiply-quantized vortex states. From Eq. (0J 
it is also clear that for harmonic trapping, A = 0, O 
cannot exceed u>, since the atoms will fly apart. For 
any other potential that grows more rapidly than p 2 at 
large distances, the effective potential is bounded, inde- 
pendently of Q. Finally, when both era and il/u> are 
sufficiently large, there is a third, mixed phase consist- 
ing of a multiply-quantized vortex state at the center of 
the cloud, surrounded by singly-quantized vortices. This 
configuration in a sense compromises between the single- 
particle energy and the interaction energy. 

Turning to the case of effective attractive interactions, 
a < 0, for small enough <r\a\, the phase diagram is, 
crudely speaking, symmetric with respect to era. The 
lowest state of the gas is determined by the single-particle 
part of the Hamiltonian, and thus the corresponding 
state still consists of multiply quantized vortex states. 
For higher values of cr\a\, these states arc again unstable 
against the formation of a combination of multiply quan- 
tized and singly quantized vortices, as in the mixed phase 
described for positive aa. The only difference is that the 
single-particle density distribution resembles that of a lo- 
calized "blob" rota ting around the minimum of the effec- 
tive potential V e g [16(. For even higher values of cr\a\, 
the phase that minimizes the energy is that of the cen- 
ter of mass excitation first seen in harmonically trapped 
atoms [22, [24| . In this phase the angular momentum is 
carried by the center of mass. The single-particle density 
distribution resembles that of the "mixed" phase (i.e., a 
localized blob) , although the structure of the many-body 
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FIG. 1: Schematic phase diagram which shows the various 
phases as the rotational frequency and the coupling are varied. 



wavefunction is very different. Fin ally, for sufficiently 
large values of cr|a| the gas collapses [2^, |2^j ■ Figure 
I shows a schematic phase diagram. 



IV. NUMERICAL DIAGONALIZATION OF THE 
HAMILTONIAN 

As mentioned earlier, in the present study we perform 
numerical diagonalization of the Hamiltonian for a given 
atom number ./V and angular momentum Lh. To achieve 
this, we first solve the single-particle eigenvalue problem 
numerically, 
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using the Numerov method. The eigenfunctions and 
eigenenergies are characterized by two quantum numbers, 
the number of radial nodes n r , and the number of quanta 
of angular momentum m. Making the simplified assump- 
tion of weak interactions, er|a| <C 1, we perform our cal- 
culation within the approximation of the lowest Landau 
level, n r = and m > 0. States with n r ^ or nega- 
tive m are higher in energy by a term of order Tiuj. This 
assumption allows us to consider relatively large values 
of N and L. One should remember that while for re- 
pulsive interactions era can have any value, for attractive 
interactions, a\a\ < 1, as otherwise the system collapses. 
Therefore, the lowest Landau level approximation is more 
restrictive for repulsive interactions, yet it gives a rather 
good description of the various phases, even in this case. 

Having found the single-particle eigenstates, we then 
set up the Fock states, which are the eigenstates of the 
number operator N and of the total angular momentum 



3 



aa=0.04, £2/(0=1 .005 



0.09 



ro 0.06 



0.03 



■ \(0,2,4 


y\o,3,6j\ 

m=2\ / \ 
Vnn=3\ 


m=4\ ■ 


X - 

m=1 


x^-X > 


H 



1.01 1.02 1.03 

Q/co 



CO 

D 



-0.002 



-0.004 



-0.006 



m=0 


jti=1 /\m=2Am=3/\m=4j 




\ xy \ x/\ x/\ x/ 




\X Mixed 
cm. \ 



1.01 1.02 

£2/co 



1.03 



FIG. 2: The phase diagram calculated within the mean- field 
approximation for A = 0.005 from Refs. |H H5 H^. The 
"cm." phase is the center of mass phase. The stars and the 
three dashed horizontal lines represent some of the values of 
a a and Sl/u> that have been examined and analyzed in the 
present study. 



L, respectively. As a final step we set up the correspond- 
ing Hamiltonian matrix between the various Fock states, 
which we diagonalize numerically. The diagonal matrix 
elements consist both of the single-particle part, as well 
as the interaction. The only non-zero off-diagonal ma- 
trix elements result from the interaction. The output of 
the Hamiltonian gives us not only the lowest state, but 
the whole energy spectrum. Since our calculation is per- 
formed for fixed angular momentum, in order to work 
with a fixed fl we consider the energy in the rotating 
frame, which is given by the usual transformation 




FIG. 3: The occupancy of the single-particle states of angular 
momentum m for A = 0.005, aa = 0.04, Q/u) = 1.005, and 
for N = 5 (higher), TV = 10 (middle), and N = 20 (lower). 
These graphs demonstrate that the occupancy of all states 
with m 7^ 1 is, to leading order, 1/N. The occupancy of 
the m = 1 state is unity minus corrections of order 1/N to 
leading order. In the mean-field approximation this phase 
corresponds to a single vortex state. 



RESULTS - REPULSIVE INTERACTIONS 



H' = H — LQ. (6) 

Having found the many-body eigenstates and eigenval- 
ues for some range of L, then for a fixed O we identify 
the eigenstate with the lowest eigenenergy according to 
Eq. 10. The many-body eigenfunction that corresponds 
to the lowest eigenenergy - which is expressed in terms 
of the Fock states - is then analyzed in terms of the oc- 
cupancy of the single particle states. The calculation of 
the occupancies is a trivial operation and, remarkably, 
each of the phases described in the previous section is 
characterized by a distinct distribution. This fact pro- 
vides indisputable evidence for the phases we expect to 
get from the arguments of the previous section and from 
the predictions of mean-field. 



We turn now to the analysis of our results, starting 
with the case of repulsive interactions. As we argued 
earlier, for sufficiently weak interactions we expect the 
phase of multiple quantization to have the lowest energy. 
For convenience we compare our results with those of the 
phase diagram of Ref. |lj|. We thus choose A = 0.005 
|28|, and evaluate the occupancy of the single particle 
states for era = 0.04, fi/w = 1.005, for N = 5, 10, and 
20 atoms, as shown in Fig. 3. This graph demonstrates 
clearly that as N increases, the occupancy of the m = 1 
state approaches unity, while the occupancy of all other 
states tends to zero (scaling as 1/N to leading order). 
In the mean-field description the order parameter ip is 
simply 

ip = 0o,i, (7) 
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FIG. 4: The occupancy of the single-particle states of angular 
momentum m for A = 0.005, aa = 0.04, N = 10 atoms, and 
Vt/tu = 1.005 (top left), Q,/lo = 1.014 (top right), M/cu = 1.020 
(bottom left), Q/cu — 1.026 (bottom right). The angular mo- 
mentum of the gas changes discontinuously between succes- 
sive values of m as SJ increases. These phases correspond to 
multiply-quantized vortex states. 



FIG. 6: The occupancy of the single-particle states of angu- 
lar momentum m for A = 0.005, aa = 0.1, fl/co = 1.0025 
and N = 10 (higher), N = 20 (lower). Here L/N = 1.8. 
These graphs demonstrate that the dominant states are the 
ones with m = 0, 2, and 4, in agreement with the mean-field 
calculations of Refs. |14Ul5|l . In the mean-field approximation 
this corresponds to a "mixed" phase discussed in the text. 



N=5, oa=0.04 




1.02 1.03 1.04 

FIG. 5: The total angular momentum of the gas in its lowest- 
energy state as a function of Q/uj, for fixed A = 0.005, 
aa = 0.04, and N — 5. The stars denote the result of a 
perturbative, mean- field calculation from Ref. |T^I . 



which represents a single vortex state located at the cen- 
ter of the cloud. 

In Fig. 4 we calculate the occupancies of the single par- 
ticle states for A = 0.005, aa — 0.04 and for four values 
of fi/w = 1.005, 1.014, 1.020, and 1.026. This graph con- 
firms that indeed the system (in its lowest state) under- 



goes discontinuous transitions between states of multiple 
quantization of successive values of m. Figure 5 shows the 
corresponding angular momentum per particle I = L/N 
in the lowest state of the gas versus Cl/u>, for a fixed value 
of aa = 0.04. This graph consists of discontinuous jumps 
in I as fl increases. The asterisks in this graph show the 
result of mean field |l5j . which compares well with our 
result for small f2. For larger f2 the agreement becomes 
worse. This is because the contribution of the anhar- 
monic term to the energy was calculated perturbatively 
in A in Ref. ^5|. This correction increases quadratically 
with m, and for larger values of this quantum number, 
the perturbative result deviates from the exact. 

For repulsive interactions, as aa increases, mean-field 
theory predicts that the system undergoes a continuous 
transition from some multiply-quantized vortex state of 
angular momentum mo to a "mixed" state which, to lead- 
ing order, is a linear superposition of three states with an- 
gular momentum mi, mo, and ma, with mi + m2 = 2mo- 
Thus the order parameter has the form 



+ c 



mo ^0,7710 



(8) 



where the coefficients are complex numbers. We have 
confirmed that our method gives a similar transition, as 
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FIG. 7: The occupancy of the single-particle states of angular 
momentum m for A = 0.005, aa — —0.001, TV = 5 atoms, and 
Q/lu = 1.013 (top left), Q,/lo = 1.018 (top right), fi/u = 1.023 
(bottom left), Q/ui = 1.028 (bottom right). In the mean-field 
approximation this phase corresponds to multiply-quantized 
vortex states. 



shown in Fig. 6. In this figure we plot the occupancy of 
the single-particle states for aa = 0.1, Vl/uj = 1.0025, 
and A = 0.005, for TV = 10 and TV = 20. From these 
two graphs we see that the states with m\ = 0, mo = 2, 
777-2 = 4 are the dominant ones for large TV. This is again 
in agreement with the prediction of mean-field. Finally, 
the angular momentum per particle I = L/N is 1.8 in 
this state. 



VI. RESULTS - ATTRACTIVE INTERACTIONS 

For an effective attractive interaction between the 
atoms the corresponding phase diagram is even richer. 
We have obtained evidence for all the phases that are ex- 
pected, namely the multiply-quantized vortex states, the 
mixed phase that consists of singly and multiply quan- 
tized vortices, the one involving center of mass excitation, 
and finally the unstable one. 

Figure 7 shows the occupancies of the single-particle 
states in the lowest state of the gas for a fixed A = 
0.005 and aa = —0.001, and for four values of Q/u) = 
1.013,1.018,1.023, and 1.028. As in Fig. 4, as ft in- 
creases, we see that states of multiple quantization (po,m. 
with successive values of m become macroscopically oc- 



FIG. 8: The angular momentum per particle in the state 
of lowest energy as a function of Q/uj, for A = 0.005, TV = 
5, and aa — —0.001 (higher), and aa = —0.002 (lower). 
The plateaus correspond to multiply-quantized vortex states, 
which become more narrow as a\a\ increases (as also shown in 
Fig. 2), in agreement with the mean-field approximation |lql . 
The small wiggles in the curves between the plateaus are an 
artifact of the discreetness of L that we have considered. 



cupied. In Fig. 8 we plot the angular momentum per par- 
ticle I = L/N versus fl/u for aa = -0.001 and -0.002. 
Here there is a difference as compared to repulsive in- 
teractions, which is also consistent with the prediction 
of mean-field theory. In order for the gas to get from 
some state of multiple quantization mo to mo + 1 (corre- 
sponding to the plateaus in Fig. 8), it has to go through 
some "mixed state" , and thus these transitions are con- 
tinuous. Furthermore, the width (in Q/uj) of the plateaus 
decreases with increasing a\a\. This effect is clearly seen 
between the two graphs in Fig. 8. 

To get evidence for the mixed phase, we increase a\a\, 
aa = —0.004, and plot the occupancies of the single parti- 
cle states in Fig. 9 for a fixed Q./u> — 1.013, and TV = 5, 10, 
and 20. Here L/N is 1.0. The dominant single-particle 
states are the ones with m = 0, 1, and 2, with a very 
small admixture of m = 3. Again, this is in agreement 
with mean-field theory, which predicts that close to the 
phase boundary, any multiply-quantized vortex state of 
strength mo is unstable against a state with three com- 
ponents of angular momentum mo — 1, mo, and mo + 1, 

i> — C mo _i0o,m o -1 + c m o 4'0,m o + Coto+1^0,to o +1- (9) 
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FIG. 9: The occupancy of the single-particle states of angular 
momentum m for A = 0.005, aa = -0.004, Q/lu = 1.013, 
and for N = 5 (higher), N = 10 (middle), and TV = 20 
(lower). Here L/N = 1.0. These graphs demonstrate that 
the dominant states are the ones with m = 0, 1, and 2, with a 
small admixture of the state with m = 3, in agreement with 
Ref. fl^l . In the mean-field approximation, this corresponds 
to the "mixed" phase discussed in the text. 



In such a linear superposition of single-particle states, the 
corresponding single-particle density distribution may 
(and it actually does) look very much like that of a local- 
ized blob |16j . although this is a phase involving vortex 
excitation. 

The phase of the center of mass excitation expected in 
this problem 0] was also clearly seen in our calculation 
for even more negative values of era. This phase was first 
discovered in harmonic traps [23L |24| . Remarkably, the 
many-body wavefunction can be written analytically and 
has a very different structure as compared to the product 
form assumed within the mean-field approximation. This 
fact makes the validity of the mean-field approximation 
questionable for attractive interactions. 

The occupancy of the single-particle states can also 
be expressed analytically. In Ref. .23] it has been shown 
that the occupancy \c rn \ 2 of a single particle state with 
angular momentum m in a many-body state with center 
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FIG. 10: The crosses show the occupancy of the single-particle 
states of angular momentum m for A = 0.005, aa = —0.05, 
TV = 5 atoms, and Q/lu = 1.012 (top left), Q/u = 1.014 
(top right), Q/uj = 1.016 (middle left), Q/w = 1.018 (middle 
right), il/u = 1.022 (bottom left), and fi/w = 1.026 (bot- 
tom right). This phase has a very large overlap with that of 
the center of mass excitation. The circles show the results 
calculated from Eq. 1101 . 

of mass excitation of L units of angular momentum and 
N atoms is 

\c m (L,N)f = !^~ 1)L "7 ! , . (10) 
N^(L — my. m\ 

The crosses in Fig. 10 show the occupancy of the single- 
particle states calculated numerically for aa — —0.05, 
A = 0.005, N = 10, and six values of fl/u> — 
1.012,1.014,1.016,1.018,1.022, and 1.026. The circles 
represent |c m | 2 of Eq. (|10fl . We confirmed that the differ- 
ence between the occupancies calculated within our study 
and the prediction of Eq. ((TU1) decreases with increasing 
N and increasing a\a\ (as long as the gas is stable). These 
data again provide clear evidence for the phase of center 
of mass excitation. 

It should also be mentioned that the phase boundary 
between the mixed phase and the one of center of mass 
excitation lies much lower than the one calculated varia- 
tionally in Ref. 0] . The reason is that the phase bound- 
ary of Ref. 0| (for these two specific phases) is just an 
upper bound, and the present calculation is consistent 
with this fact. 

Finally, the unstable phase [25L l2rl [27J also shows up 
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in our model. This is seen in the present calculation via 
the sign of the chemical potential which becomes nega- 
tive for a value of a a — 1. This result is inconsistent 
with the assumption of weak interactions, er|a| <C 1, as 
the approximation of lowest Landau level breaks down. 
To get a quantitatively accurate description of this insta- 
bility one would have to include higher Landau levels. 



weak interatomic interactions, our results are exact. We 
got evidence for all the phases that are expected, namely 
vortex excitation of multiple quantization, single quan- 
tization, and mixed; center of mass motion, and finally 
the unstable phase. The location of the phase boundaries 
were consistent with those calculated within the mean 
field approximation. 



VII. SUMMARY 

To summarize, in the present study we considered 
bosonic atoms that rotate in an anharmonic trap. We 
attacked this problem with the method of numerical di- 
agonalization of the Hamiltonian, which goes beyond the 
mean-field approximation but also forced us to consider 
small numbers of atoms. Apart from the assumption of 
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